#include <cassert>
#include <utility>
#include <cstdlib>
#include "dft.h"

using namespace std;

bool IsPowerOfTwo(int n) {
  return (n & (n - 1)) == 0;
}


int NumBits(int n) {
  assert(n > 0);
  int num_bits = 0;
  for (n -= 1; n; n >>= 1)
    num_bits++;
  return num_bits;
}


int ReverseBits(int n, int num_bits) {
  int r = 0;
  for (int i = 0; i < num_bits; ++i) {
    r <<= 1;
    r |= n & 1;
    n >>= 1;
  }
  return r;
}


void ReorderElements(complex<float>* x, int n) {
  assert(IsPowerOfTwo(n));
  int num_bits = NumBits(n);
  for (int i = 0; i < n; ++i) {
    int j = ReverseBits(i, num_bits);
    if (i < j)
      swap(x[i], x[j]);
  }
}


void DftGold(const complex<float>* x, int n, complex<float>* y) {
  for (int i = 0; i < n; ++i) {
    float ang = -2.f * static_cast<float>(M_PI) / n * i;
    y[i] = 0;
    for (int j = 0; j < n; ++j)
      y[i] += x[j] * exp(complex<float>(0, ang * j));
  }
}


void DftInplace(complex<float>* x, int n) {
  assert(IsPowerOfTwo(n));
  ReorderElements(x, n);
  for (int size = 1; size < n; size <<= 1) {
    float ang = -static_cast<float>(M_PI) / size;
    for (int i = 0; i < n; i += 2 * size) {
      for (int j = 0; j < size; ++j) {
        complex<float> a = x[i + j];
        complex<float> b = x[i + j + size];
        float ang_ = ang * j;
        b *= complex<float>(cos(ang_), sin(ang_));
        x[i + j] = a + b;
        x[i + j + size] = a - b;
      }
    }
  }
}


void FillRandom(float low, float high, int n, complex<float>* x) {
  for (int i = 0; i < n; ++i)
    x[i] = complex<float>(
        low + rand() / static_cast<float>(RAND_MAX) * (high - low),
        low + rand() / static_cast<float>(RAND_MAX) * (high - low));
}


float RelativeError(const complex<float>* gold, const complex<float>* mine, int n) {
  float err = 0.f;
  for (int i = 0; i < n; ++i)
    err = max(err, abs(gold[i] - mine[i]) / abs(gold[i]));
  return err;
}
